#!/usr/bin/python3
# -*- coding: iso-8859-1 -*-

"""
Copyright Pascal KREZEL Janvier 2012

pascal.krezel@gmail.com

TUBA est un programme informatique servant à modéliser rapidement
des tuyauteries à l'aide de SaloméMéca via un langage de commandes nommé PKL.

Ce logiciel est régi par la licence CeCILL soumise au droit français et
respectant les principes de diffusion des logiciels libres. Vous pouvez
utiliser, modifier et/ou redistribuer ce programme sous les conditions
de la licence CeCILL telle que diffusée par le CEA, le CNRS et l'INRIA 
sur le site "http://www.cecill.info".

En contrepartie de l'accessibilité au code source et des droits de copie,
de modification et de redistribution accordés par cette licence, il n'est
offert aux utilisateurs qu'une garantie limitée.  Pour les mêmes raisons,
seule une responsabilité restreinte pèse sur l'auteur du programme,  le
titulaire des droits patrimoniaux et les concédants successifs.

A cet égard  l'attention de l'utilisateur est attirée sur les risques
associés au chargement,  à l'utilisation,  à la modification et/ou au
développement et à la reproduction du logiciel par l'utilisateur étant 
donné sa spécificité de logiciel libre, qui peut le rendre complexe à 
manipuler et qui le réserve donc à des développeurs et des professionnels
avertis possédant  des  connaissances  informatiques approfondies.  Les
utilisateurs sont donc invités à charger  et  tester  l'adéquation  du
logiciel à leurs besoins dans des conditions permettant d'assurer la
sécurité de leurs systèmes et ou de leurs données et, plus généralement, 
à l'utiliser et l'exploiter dans les mêmes conditions de sécurité. 

Le fait que vous puissiez accéder à cet en-tête signifie que vous avez 
pris connaissance de la licence CeCILL, et que vous en avez accepté les
termes.


    
"""
from math import *

eps=1e-6
LesMath=[",","*","/","-",":"]   # Liste des co-caractères permettant le changement de variables



#---------------------------------------------------------------------------------------------------------------
def AddPVect(P,V) :
    return AddVect([[1,P],[1,V]]) 
    
    
#---------------------------------------------------------------------------------------------------------------
def AddSimple(V1,V2) :
    R=[u+v for u,v in zip(V1,V2)]
    return R


#---------------------------------------------------------------------------------------------------------------
def AddVect(L):
    #Somme des ki*Vi LV=[[k1,V1],[k2,V2],...]
    r=[0.0,0.0,0.0]
    for i in range(3):
        for t in L:
            k=t[0]
            V=t[1]
            r[i] += k*V[i]
    return r

#---------------------------------------------------------------------------------------------------------------
def AfficherErreur(L):
    for i in range(3):
        print
        
    print(" ERREUR:")
    
    for l in L:
        print("  ",l )
           
    for i in range(3):
        print
        
    exit(1)

#---------------------------------------------------------------------------------------------------------------
def AfficherInfo(n,T):
    if n==1 :
        print(T,)
    elif n==2 : 
        print("   ",T,)
    elif n==3:
        print("        ",T,)
    elif n==4:
        print("             ",T,)
    else :
        AfficherErreur(["Ce niveau de doc n'existe pas!"])
        
        
#---------------------------------------------------------------------------------------------------------------
# Angle formé par V2 par rapport à V1 autour de Va
def AfficherListe(L):
    print
    print
    for v in L:
        print(v)


#---------------------------------------------------------------------------------------------------------------

def AfficherSousTitre(T):
    print( 10*" "+ T)



#---------------------------------------------------------------------------------------------------------------
def AfficherTitre(n,T):
    if n==1 :
        print
        print
        print (40 * "-")
        print (40* "-")
        print(T,)
    elif n==2 : 
        print
        print
        print("   -------------------")
        print("   ",T,)
    elif n==3:
        print
        print
        print("        --------------")
        print("        ",T,)
    elif n==4:
        print
        print
        print("             ---------")
        print("             ",T,)
    else :
        AfficherErreur(["Ce niveau de doc n'existe pas!"])


#---------------------------------------------------------------------------------------------------------------
def AngleVect(V1,V2,Va):
    V1a = VectProjPlan(V1,Va)
    V2a = VectProjPlan(V2,Va)      
    v1a = VectNorm(V1a)
    nv1a = NormVect(v1a)
    v2a = VectNorm(V2a)
    nv2a = NormVect(v2a)
    if nv1a < eps or nv2a < eps:
        angle = 0
    else:
        ca = ProdScal(v1a,v2a)
        if ca>1 : ca=1
        if ca<-1 : ca=-1
        angle = acos(ca)*180/pi
        
        V3 = ProdVect(V1,V2)
        sens = ProdScal(V3,Va)
        
        if sens <  0: angle *=-1 
    
    return angle


#---------------------------------------------------------------------------------------------------------------
def AnnulerChaines(L,s):
    for v in L :
        s=s.replace(v,"")
    return s



#---------------------------------------------------------------------------------------------------------------
# Construction des vecteurs orthogonaux  vd2x,vd3x,vd4x à partir de Vd1x et Vd2x
def BaseDih(Vd1x,Vd2x):
     
    vd1x = VectNorm(Vd1x)
    vd2x = VectNorm(Vd2x)  
 
    vd3x = ProdVect(vd2x,vd1x)
    vd4x = ProdVect(vd2x,vd3x)
    
    return [vd2x,vd3x,vd4x]  

    


#---------------------------------------------------------------------------------------------------------------
def CalculerVdxi(Vd1x,Vd2x,ac,ad):

    acr =ac * pi / 180
    adr = ad * pi / 180

    v1 = VectNorm(Vd1x)
    v2 = VectNorm(Vd2x)
    v3 = ProdVect(v1,v2)
    v4 = ProdVect(v2,v3)
    #
    # Utilisation du premier angle d, dihédrique
    v5= AddVect([[-cos(adr),v4],[-sin(adr),v3]])
    #Les Vecteur tangents sont déphasés de + pi/2
    vt1 =  AddVect([[-cos(acr/2+pi/2) ,v5],[sin(acr/2+pi/2)    ,v2]])
    Vd3x =  AddVect([[-cos(acr+pi/2)   ,v5],[sin(acr+pi/2)      ,v2]])
    # Définition des nouveaux Vecteurs dihédriques courants
    acMax=180-eps
    if ac > acMax:
        Vd1x = vt1
        Vd2x = Vd3x
    else:
        Vd1x=Vd2x
        Vd2x=Vd3x

    return ([Vd1x,Vd2x])

#---------------------------------------------------------------------------------------------------------------
def CodeAppend(Code,L):  
    for v in L:
        Code.append(v)
    return 0

 
#---------------------------------------------------------------------------------------------------------------
def CodeInsert(Code,iCode,LTex):
    L=LTex[:]
    L.reverse()
    for x in L :
        Code.insert(iCode,x)

#---------------------------------------------------------------------------------------------------------------
def CodeInsertI(Code,iCode,n,LTex):
    L=LTex[:]
    L.reverse()
    for x in L:
        Code.insert(iCode," "*n*8+x)

 
#---------------------------------------------------------------------------------------------------------------
def ConvFloatTest(A):
    C=[]
    try:
        X=[float(x) for x in A]

    except:
        AfficherErreur("ERREUR : Présence de valeur non numérique après un égal")
    return 
    
#---------------------------------------------------------------------------------------------------------------
def ConvNFloat(A,n):
    C=[]
    nA = len(A)
    try:
        for i in range(nA):
            C.append(float(A[i]))
        for i in range(nA,n):
            C.append(0)
    except:
        print("")
        print("ConvNFloat : ERREUR")
        print(A)
        print("Ne peut pas être converti")
    return C

#---------------------------------------------------------------------------------------------------------------
def ConvNFloat2(A,n):
    C=[]
    nA = len(A)
    try:
        for i in range(nA):
            C.append(float(A[i]))
        for i in range(nA,n):
            C.append(0)
    except:
        print("")
        print("ConvNFloat : ERREUR")
        print(A)
        print("Ne peut pas être converti")
    return C

#---------------------------------------------------------------------------------------------------------------
def ConvNInt(A,n):
    C=[]
    nA = len(A)
    for i in range(nA):
        C.append(int(A[i]))
    for i in range(nA,n):
        C.append(0)
    return C

#---------------------------------------------------------------------------------------------------------------
def ConvNIntEval(A,n):
    C=[]
    nA = len(A)
    for i in range(nA):
        C.append(int(eval(A[i])))
    for i in range(nA,n):
        C.append(0)
    return C

#---------------------------------------------------------------------------------------------------------------
def ConvNString(nom,n):
    C=[]
    A=nom.split(",")
    nA = len(A)
    for i in range(nA):
        C.append(A[i])
    for i in range(nA,n):
        C.append("")
    return C


#---------------------------------------------------------------------------------------------------------------
def CouperListeA(L,s):
    L2=[]
    for t in L:
        if t!=s:
            L2.append(t)
        else:
            break
    return L2



#---------------------------------------------------------------------------------------------------------------
def EvalInString(S):
# Remplace les parties entre parenthèses ou entre "," par une évaluation
    S2=""
    ieval=0
    Seval=""
    
    for s in S:
        
        if ieval == 1:
            if s == ")":
                ieval=0
                u=eval(Seval)
                Seval=""
                S2 += str(u)
                continue
            else:
                Seval += s
                continue
        
        if s == "(":
            ieval = 1
        else:
            S2 += s
            
    
    
    return(S2)


#---------------------------------------------------------------------------------------------------------------
def Fractionner(v):
 
    if "#" in v:
        V=v.split("#")
        v=V[0]

    if ">" in v:
        V=v.split(">")
        v=V[0].strip()
        lien=V[1].strip()
    else:
        lien=""
        
    if ":" in v:
        V=v.split(":")
        v=V[0]
        nom=V[1].strip()
    else:
        nom=""
    
    if "=" in v:
        V=v.split("=")
        Cd=V[0].strip()
        A=V[1]      
        A = A.split(",")  
    else:
        Cd=v.strip() 
        A=[]
              
    return([Cd,A,nom,lien])


#---------------------------------------------------------------------------------------------------------------
def Fractionner2(v):
    
    if "#" in v:
        V=v.split("#")
        v=V[0]

    if ">" in v:
        V=v.split(">")
        v=V[0].strip()
        lien=V[1].strip()
    else:
        lien=""
        
    if ":" in v:
        V=v.split(":")
        v=V[0]
        nom=V[1].strip()
    else:
        nom=""
    
    if "=" in v:
        V=v.split("=")
        Cd=V[0].strip()
        A=V[1]        
    else:
        Cd=v.strip() 
        A=""
              
    return([Cd,A,nom,lien])


#---------------------------------------------------------------------------------------------------------------
def Indent(i):
    return " "*i


#---------------------------------------------------------------------------------------------------------------
def InsererList(code,i,L):
    V=RenverserListe(L)
    for x in V:
        code.insert(i,x)


#---------------------------------------------------------------------------------------------------------------
def InterPlanDroite(VPA1,VPA2,n1,n2,h):
    """
    Point d'intersection entre un plan défini par VPA1 et n1 et une droite définie par VPA2 et n2
    """ 
    [x1,y1,z1]=VPA1
    [a1,b1,c1]=n1
    [x2,y2,z2]=VPA2
    [a2,b2,c2]=n2
    
    x=((c1*c2+b1*b2+a1*a2)*x1+a1*c1*c2*h-a2*c1*c1*h+a1*b1*b2*h-a2*b1*b1*h)/(c1*c2+b1*b2+a1*a2)  
    
    y=((c1*c2+b1*b2+a1*a2)*y1+b1*c1*c2*h-b2*c1*c1*h-a1*a1*b2*h+a1*a2*b1*h)/(c1*c2+b1*b2+a1*a2)
    
    z=((c1*c2+b1*b2+a1*a2)*z1+c1*(b1*b2*h+a1*a2*h)-b1*b1*c2*h-a1*a1*c2*h)/(c1*c2+b1*b2+a1*a2)
    
    return([x,y,z])
    

#---------------------------------------------------------------------------------------------------------------
def NormVect(V) :
    Vn=VectNum(V)
    s = 0
    for x in Vn:
        s += pow(x,2)
    L = pow(s,0.5)
    return L

#---------------------------------------------------------------------------------------------------------------
def PrintListe(L):
    AfficherListe(L)

#---------------------------------------------------------------------------------------------------------------
def PrintListe(L):

    for v in L:
        print(v)


#---------------------------------------------------------------------------------------------------------------
def ProdScal(V1,V2) :
    s = 0
    for u,v in zip(V1,V2):
        s += u*v
    return s

#---------------------------------------------------------------------------------------------------------------

def ProdSimple(V1,V2) :
    R=[u*v for u,v in zip(V1,V2)]
    return R

#---------------------------------------------------------------------------------------------------------------


def ProdVect(V1,V2) :
#Produit vectoriel entre 2 vecteurs
    V3 = RotLeft( SousSimple( ProdSimple(V1,RotLeft(V2)),ProdSimple(V2,RotLeft(V1))) )
    return V3

#---------------------------------------------------------------------------------------------------------------
def Reconstruire(Cd,A,nom,lien):
    tex=Cd
    if  A != []:tex += "=" + VectString2(A)
    if nom != "":tex += ":" + nom
    if lien != "": tex += ">"+lien
    return tex

#---------------------------------------------------------------------------------------------------------------
def Reconstruire2(Cd,A,nom,lien):

    tex=Cd
    if  A != "":tex += "=" + A
    if nom != "":tex += ":" + nom
    if lien != "": tex += ">"+lien
    return tex

#---------------------------------------------------------------------------------------------------------------
def RenverserListe(L):
    R=[]
    n=len(L)
    
    for i in range(n):
        R.append(L[n-1-i])
    
    return R


#---------------------------------------------------------------------------------------------------------------
def ReplaceInListe(A,var,val):
#Remplace la lettre i par valeur x dans la chaîne de caractères S
    A2=[]
    for s in A:
        if s == var:
            A2.append(val)
        else:
            A2.append(s)
    return A2


#---------------------------------------------------------------------------------------------------------------
    """
    Remplace la lettre i par valeur x dans la chaîne de caractères S 
    si celle-ci correspond à une variable de calcul
    """
    
def ReplaceInString(S,var,val):

    S2= S[:]
    S3=""
    sval=str(val)
    nvar=len(var.strip())
    
    while 1:
        rep=0
        nS2 = len(S2)
        ivar=S2.find(var)

        if ivar != -1:
            if ivar == 0:
                if nS2 == nvar:
                    S3+=sval
                    break
                else:
                    c2=S2[nvar]
                    if c2 in LesMath:
                        S3+=sval+c2
                        S2=S2[nvar+1:]
                    else:
                        S3+=S2[:ivar+nvar]
                        S2=S2[ivar+nvar:]
                        
            else:
                c1=S2[ivar-1]
                if c1 in LesMath:
                    if ivar+nvar == nS2:
                        S3+=S2[:ivar]+sval
                        break
                    else:
                        c2=S2[ivar+nvar]
                        if c2 in LesMath:
                            S3+=S2[:ivar]+sval+c2
                            S2=S2[ivar+nvar+1:]  
                else:
                    S3+=S2
                    break          
        else:
            S3+=S2
            break 
    
    return S3



#---------------------------------------------------------------------------------------------------------------
def ReplaceInString2(S,var,val):

    S2=S[:]
    S3=""
    sval=str(val)
    sval=sval.strip()
    nvar=len(var)
    
    ivar=S2.find(var)
    while ivar!=-1:
        S3+=S2[:ivar]+sval
        S2=S2[ivar+nvar:]
        ivar=S2.find(var)
    S3 += S2
    return S3

#---------------------------------------------------------------------------------------------------------------
def ReplaceInString3(S,var,val):


    sval=str(val)
    nvar=len(var)
    
    #Remplacement à l'intérieur des parenthèses
    iPar1 = S.find("(")
    if iPar1 != -1 :
        S2 = S
        S3 = ""
        while iPar1 != -1:        
            iPar2=S2.find(")")
            iPar3=S2[iPar1+1:].find("(")
#            if iPar3!=-1:
#                if iPar3 < iPar2: AfficherErreur(["ReplaceInString3:Pas de parenthèses imbriquées dans les noms pour le moment!"])
            s=S2[iPar1+1:iPar2]
            if var in s :   
                s = ReplaceInString2(s,var,val)
                try:
                    sval2 = eval(s)
                    s2 = str(sval2)
                except:
                    s2 = "("+s+")"
            else:
                s2 = "("+s+")"
                
            S3+=S2[:iPar1]+s2
            S2=S2[iPar2+1:]
            iPar1 = S2.find("(")

        S3+=S2

       
    else:
        S3=S
 
    #Remplacement à l'extérieur des parenthèses  
    S2=S3
    S3=""
    ivar=S2.find(var)
    while ivar!=-1:
        S3+=S2[:ivar]+sval
        S2=S2[ivar+1:]
        ivar=S2.find(var)
    S3 += S2
            

        
    return S3

#-------------------------------------------------------------------------------
def RotLeft(V) :
    return([V[1],V[2],V[0]])




#-------------------------------------------------------------------------------
# Rotation de vecteur V autour de Va d'un angle c en degrés.
def RotVect(V,Va,c):

    V = VectNum(V)  
    Va = VectNum(Va)  
    va = VectNorm(Va)  
    
    b = ProdScal(V,va)
    Ua = VectRatio(va,b)
    U1 = AddVect([[1,V],[-1,Ua]])
    U2 = ProdVect(va,U1)
    
    cr = c*pi/180
    rV= AddVect([[1,Ua],[cos(cr),U1],[sin(cr),U2]])
    
    return rV

    
#-------------------------------------------------------------------------------
def SousSimple(V1,V2) :
    R=[u-v for u,v in zip(V1,V2)]
    return R


#-------------------------------------------------------------------------------
def TestColineaire(V1,V2):
    zero = 1E-10
    col=0
    if abs(V2[0])>zero:
        k=V1[0]/V2[0]
        if abs(V1[1] -V2[1]*k)<zero and abs(V1[2] -V2[2]*k)<zero: col=1
    elif abs(V2[1])>zero:
        k=V1[1]/V2[1]
        if abs(V1[0] -V2[0]*k)<zero and abs(V1[2] -V2[2]*k)<zero: col=1
    elif abs(V2[2])>zero:
        k=V1[2]/V2[2]
        if abs(V1[0] -V2[0]*k)<zero and abs(V1[1] -V2[1]*k)<zero: col=1
    return col
    
 
#---------------------------------------------------------------------------------------------------------------
def TrouverIndex(L,s):
    Lindex=[]
    nL = len(L)
    L2=L[:]
    indx=0

    while 1:
        try:
            ind=L2.index(s)
            indx+=ind
            Lindex.append(indx)
            indx+=1
            if indx == nL:break
            L2=L2[ind+1:]
        except ValueError:
            break

    return Lindex

#-------------------------------------------------------------------------------------------------------------
def TrouverLigne(L,s):
    r = 0
    for i,v in enumerate(L) :
        if s in v:
            r = i
            break
    return r

#---------------------------------------------------------------------------------------------------------------
def Vect(PA,PB) :
    return AddVect([[1,PB],[-1,PA]])   
    
    
#---------------------------------------------------------------------------------------------------------------
def VectCar(L) :
    return([str(x) for x in L])


#---------------------------------------------------------------------------------------------------------------
def VectDih(L,Vd1,Vd2,ad0,bd0):

    ad = float(ad0)
    bd = float(bd0)

    Base = BaseDih(Vd1,Vd2) 
    [vd2x,vd3x,vd4x] = Base
    
    # Utilisation du premier angle
    # Rotation autour de vd4x autour de vd2x
    Vr = RotVect(vd4x,vd2x,ad)

    # Utilisation du 2ème angle
    # Rotation de Vr autour de Vd5
    vd5x = ProdVect(vd2x, Vr)
    Vr = RotVect(Vr,vd5x,bd)    

    # Calcul du vecteur
    Vr = VectRatio(Vr,L)
    
    return Vr


#---------------------------------------------------------------------------------------------------------------
def VectInt(L) :
    try:
        Ln=[int(x) for x in L]
    except ValueError:
        print ("")
        print ("VectNum : ERREUR :")
        print ( L )
        print (" Le vecteur possede des caracteres non numeriques")
        exit(1)
    return Ln

#---------------------------------------------------------------------------------------------------------------
def VectNorm(V) :
    V2=VectNum(V)
    L = NormVect(V2)
    if L<eps :
        Vn=[0,0,0]
    else:
        Vn = [x/L for x in V2]
    return Vn


#---------------------------------------------------------------------------------------------------------------
def VectNum(L) :
    try:
        Ln=[float(x) for x in L]
    except ValueError:
        print ("")
        print ("VectNum : ERREUR :")
        print ( L )
        print (" Le vecteur possede des caracteres non numeriques")
        exit(1)
    return Ln
    

#---------------------------------------------------------------------------------------------------------------
def VectOpp(V) :
    return [x*(-1) for x in V]

#---------------------------------------------------------------------------------------------------------------
# Projection du vecteur V1 suivant le vecteur V2
def VectOrtho(V):
    [a,b,c]=V
    a=float(a)
    b=float(b)
    c=float(c)
    zero=0.0001
    Vk=[]
    if abs(a)>zero:
        if abs(b)>zero:
            Vk=[1,-a/b,0]
        elif abs(c)>zero:
            Vk=[1,0,-a/c]
        else:
            Vk=[0,1,0]
    elif abs(b)>zero:
        if abs(a)>zero:
            Vk=[-b/a,1,0]
        elif abs(c)>zero:
            Vk=[0,1,-b/c]
        else:
            Vk=[1,0,0]
    elif abs(c)>zero:
        if abs(a)>zero:
            Vk=[0,-c/a,1]
        elif abs(b)>zero:
            Vk=[-c/b,0,1]
        else:
            Vk=[0,1,0]
    return Vk


#---------------------------------------------------------------------------------------------------------------
def VectProjPlan(V1,Va) :

    V1a = ProdVect(Va,V1)
    Vp = ProdVect(V1a,Va)
    vp = VectNorm(Vp)
    k1p = ProdScal(V1,vp) 
    Vr = VectRatio(vp,k1p)
    
    return Vr

#---------------------------------------------------------------------------------------------------------------
def VectProjVect(V1,V2) :

    v2 = VectNorm(V2)
    a = ProdScal(V1,v2)
    Vr = VectRatio(v2,a)
    return Vr

#---------------------------------------------------------------------------------------------------------------
# Projection du vecteur V suivant le plan défini par Va
def VectRatio(V,d0) :

        
    d = float(d0)
    V = VectNum(V)
    L = NormVect(V)
    if L<eps :
        Vr = [0,0,0]
    else:
        Vr = [d*x/L for x in V]
    return Vr


#---------------------------------------------------------------------------------------------------------------
def VectString(A):
    tex="["
    for v in A:
        tex+= str(v)+","
        
    if tex!="[":
        tex=tex[:-1]
    
    tex +="]"
    return tex


#---------------------------------------------------------------------------------------------------------------
def VectString2(A):
    tex=""
    for v in A:
        tex+= str(v)+","
    
    if tex!="":
        tex=tex[:-1]
        
    return tex


#---------------------------------------------------------------------------------------------------------------
def VerifValeur(A) :
    rep = 1
    # Vérification des arguments
    
    try:
        for v in A :
            float(v)
    except ValueError:
        print(" Erreur : un des arguments n'est pas un chiffre.")
        rep = 0
    return rep

def VerifierArg(L) :
    n = len(L)
    rep = 1
    for v in L:
        if v == "" :
            print (" Erreur ligne ", i, " : un des arguments est vide")
            error =1
            continue
        try:
            float(v)
        except:
            rep = 0
            print (" Erreur : Valeur non numérique")
    return rep



#-------------------------------------------------------------------------------
def VerifierColinea(V1,V2) :
    v1 = VectNorm(V1)
    v2 = VectNorm(V2)
    v3 = ProdVect(v1,v2)
    nv3= NormVect(v3)
    if nv3 < eps :
        return 1
    else :
        return 0
        
#-------------------------------------------------------------------------------
def VerifierOrtho(V1,V2) :
    v1 = VectNorm(V1)
    v2 = VectNorm(V2)
    a = ProdScal(v1,v2)

    if a < eps :
        return 1
    else :
        return 0

 
#---------------------------------------------------------------------------------------------------------------
